{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 一.基本思路：马尔可夫链收敛到平稳态\n",
    "在[《12_05_PGM_马尔科夫链_初探及代码实现》](https://nbviewer.jupyter.org/github/zhulei227/ML_Notes/blob/master/notebooks/12_05_PGM_%E9%A9%AC%E5%B0%94%E7%A7%91%E5%A4%AB%E9%93%BE_%E5%88%9D%E6%8E%A2%E5%8F%8A%E4%BB%A3%E7%A0%81%E5%AE%9E%E7%8E%B0.ipynb)这一节，我们首次探索了马尔可夫链，并在本末讨论并说明了它的一个重要性质：**马尔可夫链的平稳态**，而且在后面的PageRank算法中使用了它的这个性质，这部分内容其实是要求这样一个问题：  \n",
    "\n",
    "<center>**已知马尔可夫模型的状态转移概率矩阵$P$，求它的平稳态$\\pi^*$，使得$P\\pi^*=\\pi^*$**</center>  \n",
    "\n",
    "将这个问题反过来，那就是**MCMC**（马尔科夫蒙特卡洛抽样法）的主要思想咯：   \n",
    "\n",
    "<center>**已知某分布$\\pi$，求一个马尔可夫模型的状态转移概率矩阵$P^*$，使得$P^*\\pi=\\pi$**</center>    \n",
    "\n",
    "于是，我们求一个$P^*$满足上面的条件即可，但是要用于MCMC，对$P^*$的要求更加苛刻一些：要保证在$P^*$的条件下，**平稳态必须是唯一**的，换言之，给定一个概率转移矩阵，它的平稳态可能有多个，比如下面的例子"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "os.chdir('../')\n",
    "from ml_models.pgm import SimpleMarkovModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "smm=SimpleMarkovModel(status_num=3)\n",
    "smm.P=np.asarray([\n",
    "    [1,1/3,0],\n",
    "    [0,1/3,0],\n",
    "    [0,1/3,1]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>50</th>\n",
       "      <th>100</th>\n",
       "      <th>500</th>\n",
       "      <th>1000</th>\n",
       "      <th>2000</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.433333</td>\n",
       "      <td>0.549994</td>\n",
       "      <td>5.500000e-01</td>\n",
       "      <td>5.500000e-01</td>\n",
       "      <td>5.500000e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.233333</td>\n",
       "      <td>0.000012</td>\n",
       "      <td>9.750689e-25</td>\n",
       "      <td>1.358228e-48</td>\n",
       "      <td>2.635403e-96</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.333333</td>\n",
       "      <td>0.449994</td>\n",
       "      <td>4.500000e-01</td>\n",
       "      <td>4.500000e-01</td>\n",
       "      <td>4.500000e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       50        100           500           1000          2000\n",
       "0  0.433333  0.549994  5.500000e-01  5.500000e-01  5.500000e-01\n",
       "1  0.233333  0.000012  9.750689e-25  1.358228e-48  2.635403e-96\n",
       "2  0.333333  0.449994  4.500000e-01  4.500000e-01  4.500000e-01"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#定义初始状态\n",
    "init_prob=np.asarray([[0.2],[0.7],[0.1]])\n",
    "pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>50</th>\n",
       "      <th>100</th>\n",
       "      <th>500</th>\n",
       "      <th>1000</th>\n",
       "      <th>2000</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.833333</td>\n",
       "      <td>0.849999</td>\n",
       "      <td>8.500000e-01</td>\n",
       "      <td>8.500000e-01</td>\n",
       "      <td>8.500000e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.033333</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>1.392956e-25</td>\n",
       "      <td>1.940325e-49</td>\n",
       "      <td>3.764862e-97</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.133333</td>\n",
       "      <td>0.149999</td>\n",
       "      <td>1.500000e-01</td>\n",
       "      <td>1.500000e-01</td>\n",
       "      <td>1.500000e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       50        100           500           1000          2000\n",
       "0  0.833333  0.849999  8.500000e-01  8.500000e-01  8.500000e-01\n",
       "1  0.033333  0.000002  1.392956e-25  1.940325e-49  3.764862e-97\n",
       "2  0.133333  0.149999  1.500000e-01  1.500000e-01  1.500000e-01"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#定义另一组初始状态\n",
    "init_prob=np.asarray([[0.8],[0.1],[0.1]])\n",
    "pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>50</th>\n",
       "      <th>100</th>\n",
       "      <th>500</th>\n",
       "      <th>1000</th>\n",
       "      <th>2000</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.566667</td>\n",
       "      <td>0.649996</td>\n",
       "      <td>6.500000e-01</td>\n",
       "      <td>6.500000e-01</td>\n",
       "      <td>6.500000e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.166667</td>\n",
       "      <td>0.000008</td>\n",
       "      <td>6.964778e-25</td>\n",
       "      <td>9.701626e-49</td>\n",
       "      <td>1.882431e-96</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.266667</td>\n",
       "      <td>0.349996</td>\n",
       "      <td>3.500000e-01</td>\n",
       "      <td>3.500000e-01</td>\n",
       "      <td>3.500000e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       50        100           500           1000          2000\n",
       "0  0.566667  0.649996  6.500000e-01  6.500000e-01  6.500000e-01\n",
       "1  0.166667  0.000008  6.964778e-25  9.701626e-49  1.882431e-96\n",
       "2  0.266667  0.349996  3.500000e-01  3.500000e-01  3.500000e-01"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#再定义另一组初始状态\n",
    "init_prob=np.asarray([[0.4],[0.5],[0.1]])\n",
    "pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),\n",
    "              2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上面的例子可以发现，由于初始状态的不同，平稳态可能有多个，所以我们构造的$P$必须要保证平稳态能唯一的收敛到我们的目标分布$\\pi$，不然，如上面的例子，假如我们的目标分布是$\\pi=[0.65,0,0.35]$，然后我们构造了一个概率转移矩阵$P=\\begin{bmatrix}\n",
    "1 & 1/3 & 0\\\\ \n",
    "0 & 1/3 &0 \\\\ \n",
    "0 & 1/3 & 1\n",
    "\\end{bmatrix}$ ，满足$P\\pi=\\pi$，然后，我们随便定义了一个初始点$x_0=[0.8,0.1,0.1]$，在$P$上随机游走采样$\\{x_1,x_2,...,x_n\\}$，最后发现它成功收敛到了$x_n\\rightarrow [0.85,0,0.15]$，哈哈哈哈哈~~~~     \n",
    "\n",
    "保证马尔科夫模型仅收敛到唯一平稳态的**充分条件**是有的！那就是**细致平衡方程**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 二.保证平稳态唯一：细致平衡方程\n",
    "细致平衡方程的定义如下，设有马尔科夫链$X=\\{X_0,X_1,...,X_t,...\\}$，状态空间为$S$，转移概率矩阵为$P$，如果有状态分布$\\pi=(\\pi_1,\\pi_2,...)^T$，对任意时刻状态$i,j\\in S$，对任意时刻$t$均满足：   \n",
    "\n",
    "$$\n",
    "P(X_t=i\\mid X_{t-1}=j)\\pi_j=P(X_{t-1}=j\\mid X_t=i)\\pi_i,i,j=1,2,....\n",
    "$$  \n",
    "\n",
    "或者简写为：   \n",
    "\n",
    "$$\n",
    "p_{ij}\\pi_j=p_{ji}\\pi_i,i,j=1,2,....\n",
    "$$    \n",
    "\n",
    "为了更加直观，也可以写作：    \n",
    "\n",
    "$$\n",
    "p(j\\rightarrow i)\\pi_j=p(i\\rightarrow j)\\pi_i\n",
    "$$\n",
    "这便是细致平衡方程，此时的马尔科夫链称为可逆马尔科夫链，显然对$P$矩阵按列求和为1，即$\\sum_{i}p_{ij}=1$或者$\\sum_ip(j\\rightarrow i)=1$ ，接下来简单证明一下，满足细致平衡条件，则有$P\\pi=\\pi$成立：   \n",
    "\n",
    "$$\n",
    "(P\\pi)_i=\\sum_jp_{ij}\\pi_j=\\sum_jp_{ji}\\pi_i=\\pi_i\\sum_{j}p_{ji}=\\pi_i,i=1,2,...\n",
    "$$   \n",
    "\n",
    "**高维/连续状态空间**：上面只是对一维离散状态空间做的说明，而细致平衡方程对高维空间或者连续状态空间一样是成立的，只需对相应符号做修改即可，比如将$p_{ij}$修改为一个函数的形式$p(状态j,状态i)$，求和符号$\\sum$需要替换为积分符号$\\int$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 三.小结\n",
    "最后小结一下，这一节的重点引出细致平衡方程：   \n",
    "\n",
    "$$\n",
    "p_{ij}\\pi_j=p_{ji}\\pi_i\n",
    "$$\n",
    "\n",
    "因为在满足该方程约束的条件下，我们构造的$P^*$，必然有$P^*\\pi=\\pi$，且$\\pi$唯一（这里$\\pi$是我们的目标分布），这也是检验MCMC是否有效的一条重要标准，下一节将要介绍的Metropolis-Hastings算法（MH算法）便是满足该约束条件的一套算法\n",
    "\n",
    "#### 采样步骤\n",
    "\n",
    "再对MCMC的采样步骤做个整理    \n",
    "\n",
    "（1）首先在随机变量$x$的状态空间上构造一个马尔科夫链$P$，使它能平稳且唯一的收敛到我们的目标分布$p(x)$；    \n",
    "\n",
    "（2）从状态空间某一点$x_0$出发，用构造的马尔科夫链$P$进行随机游走，产生样本序列$x_0,x_1,...,x_t,....$；   \n",
    "\n",
    "（3）确定一个正整数$m,n$（$m<n$），得到样本集合$\\{x_m,x_{m+1},...,x_n\\}$，对目标指标进行计算，比如某函数$f(x)$的期望：    \n",
    "\n",
    "$$\n",
    "\\hat{E}f=\\frac{1}{n-m}\\sum_{i=m+1}^nf(x_i)\n",
    "$$\n",
    "\n",
    "##### 补充\n",
    "\n",
    "（1）上面提到的**随机游走采样**简单解释一下，它其实和[《12_06_PGM_马尔科夫链_语言模型及文本生成》](https://nbviewer.jupyter.org/github/zhulei227/ML_Notes/blob/master/notebooks/12_06_PGM_%E9%A9%AC%E5%B0%94%E7%A7%91%E5%A4%AB%E9%93%BE_%E8%AF%AD%E8%A8%80%E6%A8%A1%E5%9E%8B%E5%8F%8A%E6%96%87%E6%9C%AC%E7%94%9F%E6%88%90.ipynb)中的文本生成类似，只是每一步不会做贪婪搜索或者beam search，而是按照当前步的概率分布进行采样，例如当前样本在状态$i$，则按照$p_{ji},j=1,2,...$的分布对下一个样本进行采样；     \n",
    "\n",
    "（2）另外由于马尔科夫链的性质，当$t\\rightarrow \\infty$时，采样的样本分布才会更加近似于目标分布，所以通常我们会从所有的采样样本$\\{x_0,x_1,x_2,...,x_m,...,x_n\\}$中截取尾部的一部分样本$\\{x_m,x_{m+1},...,x_n\\}$来做相应的统计指标计算，比如计算期望，积分等...   \n",
    "\n",
    "（3）**另外另外**，有时为了方便直观理解，转移概率我有时会采用不同的写法，但请注意一下它们是等价的：$p_{ij}=p(i\\mid j)=p(j\\rightarrow i)=p(j,i)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
